This is code for Section 17.2 of Kruschke using a Student t noise distribution. Other than that change, this document is very similar to the one from section 7.1.
Make sure that the two files “HtWtData30.csv” and “HtWtData300.csv” are in your project directory.
Load necessary packages:
library(tidyverse)
library(rstan)
library(tidybayes)
library(bayesplot)
Set Stan to save compiled code.
rstan_options(auto_write = TRUE)
Set Stan to use parallel processing where possible.
options(mc.cores = parallel::detectCores())
HtWtData30 <- read_csv("HtWtData30.csv")
Rows: 30 Columns: 3── Column specification ─────────────────────────────────────────
Delimiter: ","
dbl (3): male, height, weight
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
HtWtData30
ggplot(HtWtData30, aes(y = weight, x = height)) +
geom_point()
HtWtData300 <- read_csv("HtWtData300.csv")
Rows: 300 Columns: 3── Column specification ─────────────────────────────────────────
Delimiter: ","
dbl (3): male, height, weight
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
HtWtData300
ggplot(HtWtData300, aes(y = weight, x = height)) +
geom_point()
Mean center the x variable. (Kruschke standardizes both variables using z-scores, but that makes the coefficients of the model harder to interpret.)
HtWtData30 <- HtWtData30 %>%
mutate(height_mc = height - mean(height))
HtWtData30
ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
geom_point()
HtWtData300 <- HtWtData300 %>%
mutate(height_mc = height - mean(height))
HtWtData300
ggplot(HtWtData300, aes(y = weight, x = height_mc)) +
geom_point()
We’re going to run two different examples, so we’ll make two lists.
stan_data_30 <- list(N = NROW(HtWtData30),
y = HtWtData30$weight,
x = HtWtData30$height_mc)
stan_data_300 <- list(N = NROW(HtWtData300),
y = HtWtData300$weight,
x = HtWtData300$height_mc)
data {
int<lower = 0> N;
vector<lower = 0>[N] y;
vector[N] x;
}
generated quantities {
real nu;
real beta0;
real beta1;
vector[N] mu;
real<lower = 0> sigma;
real y_sim[N];
nu = gamma_rng(2, 0.1); // default prior for df
beta0 = normal_rng(150, 25); // Intercept between 100 and 200
beta1 = normal_rng(0, 5); // -10 to 10 lbs per inch
mu = beta0 + beta1 * x; // linear model
sigma = exponential_rng(0.1); // default prior for sigma
y_sim = student_t_rng(nu, mu, sigma);
}
fit_HtWt30_robust_prior <- sampling(HtWt_robust_prior,
data = stan_data_30,
chains = 1,
algorithm = "Fixed_param",
seed = 11111,
refresh = 0)
samples_HtWt30_prior <- tidy_draws(fit_HtWt30_robust_prior)
samples_HtWt30_prior
y_sim_prior <- samples_HtWt30_prior %>%
select(starts_with("y_sim")) %>%
as.matrix()
ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
geom_point() +
geom_abline(data = samples_HtWt30_prior,
aes(intercept = beta0, slope = beta1),
alpha = 0.05)
mcmc_hist(samples_HtWt30_prior,
pars = c("nu", "sigma", "beta0", "beta1"))
ppc_intervals(y = HtWtData30$weight,
x = HtWtData30$height,
yrep = y_sim_prior)
data {
int<lower = 0> N;
vector<lower = 0>[N] y;
vector[N] x;
}
parameters {
real nu;
real beta0;
real beta1;
real<lower = 0> sigma;
}
transformed parameters {
vector[N] mu;
mu = beta0 + beta1 * x; // linear model
}
model {
nu ~ gamma(2, 0.1); // default prior for df
beta0 ~ normal(150, 25); // Intercept between 100 and 200
beta1 ~ normal(0, 5); // -10 to 10 lbs per inch
sigma ~ exponential(0.1); // default prior for sigma
y ~ student_t(nu, mu, sigma);
}
generated quantities {
real y_sim[N];
y_sim = student_t_rng(nu, mu, sigma);
}
fit_HtWt30_robust <- sampling(HtWt_robust,
data = stan_data_30,
seed = 11111,
refresh = 0)
Warning: There were 676 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
plot(fit_HtWt30_robust, plotfun = "ac",
pars = c("nu", "beta0", "beta1", "sigma"))
plot(fit_HtWt30_robust, plotfun = "trace",
pars = c("nu", "beta0", "beta1", "sigma"))
print(fit_HtWt30_robust, pars = c("nu", "beta0", "beta1", "sigma"))
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5%
nu 16.17 0.49 11.92 2.95 7.17 12.86 21.76 47.73
beta0 152.19 0.13 4.72 143.17 149.05 152.09 155.14 161.87
beta1 4.31 0.04 1.21 1.85 3.53 4.33 5.12 6.55
sigma 24.05 0.19 4.19 16.65 21.08 23.79 26.67 32.84
n_eff Rhat
nu 584 1
beta0 1279 1
beta1 1191 1
sigma 481 1
Samples were drawn using NUTS(diag_e) at Fri May 12 12:08:46 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
pairs(fit_HtWt30_robust,
pars = c("nu", "beta0", "beta1", "sigma"))
stan_dens(fit_HtWt30_robust,
pars = c("nu", "beta0", "beta1", "sigma"))
plot(fit_HtWt30_robust, pars = c("nu", "beta0", "beta1", "sigma"))
ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)
fit_HtWt300_robust <- sampling(HtWt_robust,
data = stan_data_300,
seed = 11111,
refresh = 0)
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problems
plot(fit_HtWt300_robust, plotfun = "ac",
pars = c("nu", "beta0", "beta1", "sigma"))
plot(fit_HtWt300_robust, plotfun = "trace",
pars = c("nu", "beta0", "beta1", "sigma"))
print(fit_HtWt300_robust,
pars = c("nu", "beta0", "beta1", "sigma"))
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5%
nu 5.59 0.05 1.79 3.27 4.34 5.25 6.39 10.11
beta0 156.24 0.03 1.67 153.03 155.10 156.22 157.36 159.52
beta1 4.43 0.01 0.41 3.60 4.15 4.43 4.70 5.22
sigma 23.95 0.04 1.66 20.84 22.77 23.92 25.01 27.46
n_eff Rhat
nu 1142 1
beta0 3146 1
beta1 3097 1
sigma 1656 1
Samples were drawn using NUTS(diag_e) at Fri May 12 12:18:30 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
pairs(fit_HtWt300_robust,
pars = c("nu", "beta0", "beta1", "sigma"))
plot(fit_HtWt300_robust, plotfun = "dens",
pars = c("nu", "beta0", "beta1", "sigma"))
plot(fit_HtWt300_robust,
pars = c("nu", "beta0", "beta1", "sigma"))
ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)
We’ll just do this for the example with 30 observations.
samples_HtWt30 <- tidy_draws(fit_HtWt30_robust)
samples_HtWt30
y_sim_post <- samples_HtWt30 %>%
select(starts_with("y_sim")) %>%
as.matrix()
ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
geom_point() +
geom_abline(data = samples_HtWt30,
aes(intercept = beta0, slope = beta1),
alpha = 0.01) +
geom_abline(data = samples_HtWt30,
aes(intercept = mean(beta0), slope = mean(beta1)),
color = "blue", size = 2)
Compare to standard linear regression.
ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
geom_point() +
geom_abline(data = samples_HtWt30,
aes(intercept = beta0, slope = beta1),
alpha = 0.01) +
geom_abline(data = samples_HtWt30,
aes(intercept = mean(beta0), slope = mean(beta1)),
color = "blue", size = 2) +
geom_smooth(method = lm, color = "red", size = 2)
ppc_intervals(y = HtWtData30$weight,
x = HtWtData30$height,
yrep = y_sim_post)
ppc_hist(y, y_sim_post[1:5, ])
ppc_boxplot(y, y_sim[1:5, ])
ppc_dens(y, y_sim[1:5, ])
ppc_dens_overlay(y, y_sim[1:30, ])
ppc_stat_2d(y, y_sim)